///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/wfn/symmetry.h
///
/// \brief   header file for point group symmetry
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef FORMIC_WFN_SYMMETRY_HEADER
#define FORMIC_WFN_SYMMETRY_HEADER

#include<vector>
#include<string>

#include<boost/shared_ptr.hpp>

#include<formic/exception.h>

namespace formic {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  provides the direct product table for the requested point group
  ///
  /// \param[in]      gname    the point group's name
  /// \param[out]     table    on exit, the direct product table
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  inline void group_direct_product_table(const std::string & gname, std::vector<int> & table) {

    if ( gname == "C1" ) {

      table.assign(1, 0);

    } else if ( gname == "Ci" || gname == "C2" || gname == "Cs" ) {

      // Ci ordering:  0   1
      // Ci ordering:  Ag  Au

      // C2 ordering:  0   1
      // C2 ordering:  A   B

      // C2 ordering:  0   1
      // C2 ordering:  A'  A''

      const int n = 2;
      table.resize(n*n);
      table.at(0+n*0) = 0;  table.at(0+n*1) = 1;
      table.at(1+n*0) = 1;  table.at(1+n*1) = 0;

    } else if ( gname == "C2v" || gname == "C2h" || gname == "D2" ) {

      // C2v  ordering:  0   1   2   3
      // C2v  ordering:  A1  A2  B1  B2

      // C2h  ordering:  0   1   2   3
      // C2h  ordering:  Ag  Bg  Au  Bu

      // D2   ordering:  0   1   2   3
      // D2   ordering:  A   B1  B2  B3

      const int n = 4;
      table.resize(n*n);
      table.at(0+n*0) = 0;  table.at(0+n*1) = 1;  table.at(0+n*2) = 2;  table.at(0+n*3) = 3;
      table.at(1+n*0) = 1;  table.at(1+n*1) = 0;  table.at(1+n*2) = 3;  table.at(1+n*3) = 2;
      table.at(2+n*0) = 2;  table.at(2+n*1) = 3;  table.at(2+n*2) = 0;  table.at(2+n*3) = 1;
      table.at(3+n*0) = 3;  table.at(3+n*1) = 2;  table.at(3+n*2) = 1;  table.at(3+n*3) = 0;

    } else if ( gname == "D2h" ) {

      // ordering:    0     1     2     3     4     5     6     7
      // ordering:    Ag    B1g   B2g   B3g   Au    B1u   B2u   B3u
      const int n = 8;
      std::vector<int> & t = table;
      t.resize(n*n);
      t[0+n*0] = 0;  t[0+n*1] = 1;  t[0+n*2] = 2;  t[0+n*3] = 3;  t[0+n*4] = 4;  t[0+n*5] = 5;  t[0+n*6] = 6;  t[0+n*7] = 7;
      t[1+n*0] = 1;  t[1+n*1] = 0;  t[1+n*2] = 3;  t[1+n*3] = 2;  t[1+n*4] = 5;  t[1+n*5] = 4;  t[1+n*6] = 7;  t[1+n*7] = 6;
      t[2+n*0] = 2;  t[2+n*1] = 3;  t[2+n*2] = 0;  t[2+n*3] = 1;  t[2+n*4] = 6;  t[2+n*5] = 7;  t[2+n*6] = 4;  t[2+n*7] = 5;
      t[3+n*0] = 3;  t[3+n*1] = 2;  t[3+n*2] = 1;  t[3+n*3] = 0;  t[3+n*4] = 7;  t[3+n*5] = 6;  t[3+n*6] = 5;  t[3+n*7] = 4;
      t[4+n*0] = 4;  t[4+n*1] = 5;  t[4+n*2] = 6;  t[4+n*3] = 7;  t[4+n*4] = 0;  t[4+n*5] = 1;  t[4+n*6] = 2;  t[4+n*7] = 3;
      t[5+n*0] = 5;  t[5+n*1] = 4;  t[5+n*2] = 7;  t[5+n*3] = 6;  t[5+n*4] = 1;  t[5+n*5] = 0;  t[5+n*6] = 3;  t[5+n*7] = 2;
      t[6+n*0] = 6;  t[6+n*1] = 7;  t[6+n*2] = 4;  t[6+n*3] = 5;  t[6+n*4] = 2;  t[6+n*5] = 3;  t[6+n*6] = 0;  t[6+n*7] = 1;
      t[7+n*0] = 7;  t[7+n*1] = 6;  t[7+n*2] = 5;  t[7+n*3] = 4;  t[7+n*4] = 3;  t[7+n*5] = 2;  t[7+n*6] = 1;  t[7+n*7] = 0;

    } else {
      throw formic::Exception("unrecognized point group name \"%s\" in formic::group_direct_product_table") % gname;
    }

  }

}

#endif
